/*
 * GridSpace.cpp
 *
 *  Created on: Nov 9, 2011
 *      Author: yi
 */

#include "mbuild/GridSpace.h"

GridSpace::GridSpace() {
	// initialize an empty GridSpace with size zero.
	x0=y0=z0=0;
	ncel=ncel_x=ncel_y=ncel_z=0;
	dx=dy=dz=0;
	xshift=yshift=zshift=0;
	data=NULL;
}

GridSpace::~GridSpace() {
	delete [] data;
}

void GridSpace::set_dimension(unsigned long nx, unsigned long ny, unsigned long nz) {
	this->ncel_x=nx;
//DEBUG	printf("%lu,%lu,%lu\n",ULONG_MAX,ULONG_MAX/nx,ULONG_MAX/(nx*ny));
	if(ULONG_MAX/nx <= ny)  throw TooBigForThisPlatform();
	this->ncel_y=ny;
	if(ULONG_MAX/(nx*ny) <= nz) throw TooBigForThisPlatform();
	this->ncel_z=nz;
	this->ncel=nx*ny*nz;
}



unsigned long GridSpace::get_dimension(unsigned long & nx, unsigned long & ny, unsigned long & nz) {
	nx=this->ncel_x;
	ny=this->ncel_y;
	nz=this->ncel_z;
	return this->ncel;
}



void GridSpace::set_origin(double x0, double y0, double z0) {
	/*
	 * Set the origin(the lowest corner of the box).
	 */
	this->x0=x0;
	this->y0=y0;
	this->z0=z0;
}



void GridSpace::get_origin(double & x0, double & y0, double & z0) {
	x0=this->x0;
	y0=this->y0;
	z0=this->z0;
}



void GridSpace::set_grid_size(double dx, double dy, double dz)
{
	if(dx >0 && dy >0 && dz >0) {
	this->dx=dx;
	this->dy=dy;
	this->dz=dz;
	} else {
		throw GridSizeIsNotPositive();
	}
}



void GridSpace::get_grid_size(double & dx, double & dy, double & dz) {
	dx=this->dx;
	dy=this->dy;
	dz=this->dz;
}


double* GridSpace::values(void) {
	/*
	 *  Return the pointer for the value array.
	 */
	return this->data;
}

void GridSpace::initialize(double value) {
	/*
	 *  Allocate space in memory and initialize the space with an initial value
	 */
	unsigned long offset=0;
	// If already initialized, delete the data and re-initialize.
	if (this->data != NULL) delete [] this->data;
	try {
	this->data= new double[ncel];
	} catch(bad_alloc) {
		throw NotEnoughMemory();
	}
	if (this->data == NULL) {
		throw NotEnoughMemory();
	}
	// give each grid an initial value.
	double x1=x0+xshift;
	double y1=y0+yshift;
	double z1=z0+zshift;
	double xh=x1+ncel_x*dx;
	double yh=y1+ncel_y*dy;
	double zh=z1+ncel_z*dz;
	for(int i=0;i<ncel_x;i++)
	{
	  for(int j=0;j<ncel_y;j++)
		 {
			 for(int k=0;k<ncel_z;k++)
			{
				offset=ncel_y*ncel_z*i+ncel_z*j+k;
				data[offset]=value;
			}
		}
	}
}

void GridSpace::clear(double value) {
	/*
	 *  Clear the values of the grids.
	 */
	  for(register unsigned long i=0;i<ncel;i++)
	    {
	      data[i]=value;
	    }
}




double GridSpace::getDerivatives(double x, double y, double z, double & dvdx, double & dvdy, double & dvdz,BOUNDARY bnd) {
	/*
	 * Get the first order derivatives (dvdx,dvdy,dvdz) using tri-linear interpolation.
	 */
	double ixf,iyf,izf,dx,dy,dz;
	double value;
	int ix,iy,iz;
	unsigned long index;
	if(this->dx== 0 || this->dy ==0 || this->dz==0) throw GridSpaceNotInitialized();
	ixf=(x-x0)/this->dx;
	iyf=(y-y0)/this->dy;
	izf=(z-z0)/this->dz;
	ix=floor(ixf);
	iy=floor(iyf);
	iz=floor(izf);
	dx=ixf-(double)(ix);
	dy=iyf-(double)(iy);
	dz=izf-(double)(iz);
	//If x,y,z is not inside the box
	//TODO: Improve maybe? Periodic boundary.
	switch(bnd){
	    case NO : {
	        if(ix<0 || ix>=ncel_x || iy<0 || iy>=ncel_y || iz<0 || iz>=ncel_z) throw OutOfRange();
	    }
	    case PEROIDIC : {
	        if(ix<0) ix=ncel_x-1 - (-ix) % (ncel_x-1);
	        if(ix>= ncel_x) ix=ix % (ncel_x-1);
	        if(iy<0) iy=ncel_y-1- (-iy) % (ncel_y-1);
	        if(iy>= ncel_y) iy=iy % (ncel_y-1);
	        if(iz<0) iz=ncel_z-1 - (-iz) % (ncel_z-1);
	        if(iz>= ncel_z) iz=iz % (ncel_z-1);
	    }
	}

	//Calc value

    index=ix*this->ncel_x*this->ncel_y+iy*this->ncel_z+iz;
    value=dx*dy*dz*this->data[index+ncel_y*ncel_z+ncel_z+1];
    value+=dx*dy*(1-dz)*this->data[index+ncel_y*ncel_z+ncel_z];
    value+=dx*(1-dy)*dz*this->data[index+ncel_y*ncel_z+1];
    value+=dx*(1-dy)*(1-dz)*this->data[index+ncel_y*ncel_z];
    value+=(1-dx)*dy*dz*this->data[index+ncel_z+1];
    value+=(1-dx)*dy*(1-dz)*this->data[index+ncel_z];
    value+=(1-dx)*(1-dy)*dz*this->data[index+1];
    value+=(1-dx)*(1-dx)*(1-dz)*this->data[index];

    //Calc dv/dx;
    dvdx=dy*dz*this->data[index+ncel_y*ncel_z+ncel_z+1];
    dvdx+=dy*(1-dz)*this->data[index+ncel_y*ncel_z+ncel_z];
    dvdx+=(1-dy)*dz*this->data[index+ncel_y*ncel_z+1];
    dvdx+=(1-dy)*(1-dz)*this->data[index+ncel_y*ncel_z];
    dvdx-=dy*dz*this->data[index+ncel_z+1];
    dvdx-=dy*(1-dz)*this->data[index+ncel_z];
    dvdx-=(1-dy)*dz*this->data[index+1];
    dvdx-=(1-dy)*(1-dz)*this->data[index];
    dvdx=dvdx/this->dx;
    //Calc dv/dy;
    dvdy=dx*dz*this->data[index+ncel_y*ncel_z+ncel_z+1];
    dvdy+=dx*(1-dz)*this->data[index+ncel_y*ncel_z+ncel_z];
    dvdy-=dx*dz*this->data[index+ncel_y*ncel_z+1];
    dvdy-=dx*(1-dz)*this->data[index+ncel_y*ncel_z];
    dvdy+=(1-dx)*dz*this->data[index+ncel_z+1];
    dvdy+=(1-dx)*(1-dz)*this->data[index+ncel_z];
    dvdy-=(1-dx)*dz*this->data[index+1];
    dvdy-=(1-dx)*(1-dz)*this->data[index];
    dvdy=dvdy/this->dy;
    //Calc dv/dz;
    dvdz=dx*dy*this->data[index+ncel_y*ncel_z+ncel_z+1];
    dvdz-=dx*dy*this->data[index+ncel_y*ncel_z+ncel_z];
    dvdz+=dx*(1-dy)*this->data[index+ncel_y*ncel_z+1];
    dvdz-=dx*(1-dy)*this->data[index+ncel_y*ncel_z];
    dvdz+=(1-dx)*dy*this->data[index+ncel_z+1];
    dvdz-=(1-dx)*dy*this->data[index+ncel_z];
    dvdz+=(1-dx)*(1-dy)*this->data[index+1];
    dvdz-=(1-dx)*(1-dy)*this->data[index];
    dvdz=dvdz/this->dz;
    return value;
}

double GridSpace::getValue(double x, double y, double z, BOUNDARY bnd) {
	/*
	 * Get the value of data point(x,y,z) from tri-linear inteporation.
	 */
	double ixf,iyf,izf,dx,dy,dz;
	double value;
	int ix,iy,iz;
	unsigned long index;
	if(this->dx== 0 || this->dy ==0 || this->dz==0) throw GridSpaceNotInitialized();
	ixf=(x-x0)/this->dx;
	iyf=(y-y0)/this->dy;
	izf=(z-z0)/this->dz;
	ix=floor(ixf);
	iy=floor(iyf);
	iz=floor(izf);
	dx=ixf-(double)(ix);
	dy=iyf-(double)(iy);
	dz=izf-(double)(iz);
	//If x,y,z is not inside the box
	switch(bnd){
	    case NO : {
	        if(ix<0 || ix>=ncel_x-1 || iy<0 || iy>=ncel_y-1 || iz<0 || iz>=ncel_z-1) throw OutOfRange();
	    }
	    case PEROIDIC : {
	        if(ix<0) ix=ncel_x-1 - (-ix) % (ncel_x-1);
	        if(ix>= ncel_x) ix=ix % (ncel_x-1);
	        if(iy<0) iy=ncel_y-1- (-iy) % (ncel_y-1);
	        if(iy>= ncel_y) iy=iy % (ncel_y-1);
	        if(iz<0) iz=ncel_z-1 - (-iz) % (ncel_z-1);
	        if(iz>= ncel_z) iz=iz % (ncel_z-1);
	    }
	}
	// Calc value;
    index=ix*this->ncel_x*this->ncel_y+iy*this->ncel_z+iz;
    value=dx*dy*dz*this->data[index+ncel_y*ncel_z+ncel_z+1];
    value+=dx*dy*(1-dz)*this->data[index+ncel_y*ncel_z+ncel_z];
    value+=dx*(1-dy)*dz*this->data[index+ncel_y*ncel_z+1];
    value+=dx*(1-dy)*(1-dz)*this->data[index+ncel_y*ncel_z];
    value+=(1-dx)*dy*dz*this->data[index+ncel_z+1];
    value+=(1-dx)*dy*(1-dz)*this->data[index+ncel_z];
    value+=(1-dx)*(1-dy)*dz*this->data[index+1];
    value+=(1-dx)*(1-dx)*(1-dz)*this->data[index];
    return value;
}

inline double GridSpace::value(int ix, int iy, int iz)
{
	if(ncel_x == 0 || ncel_y == 0 || ncel_z==0) throw GridSpaceNotInitialized();
	if( ix<0 || ix>=ncel_x || \
	    iy<0 || iy>=ncel_y || \
	    iz<0 || iz>=ncel_z) throw OutOfRange();
	return data[ix*ncel_y*ncel_z+iy*ncel_z+iz];
}

void GridSpace::copyto(GridSpace *target) {
	/*
	 * Copy the current data to a target Grid Space.
	 */
    target->set_dimension(ncel_x,ncel_y,ncel_z);
    target->set_grid_size(dx,dy,dz);
    target->set_origin(x0,y0,z0);
    target->initialize();
    for(unsigned long i=0;i<ncel;i++) {
       target->data[i]=this->data[i];
    }
}

bool GridSpace::eql(GridSpace & g)
{
	unsigned long nnx,nny,nnz,nn;
	double ndx,ndy,ndz,nx0,ny0,nz0;
	nn=g.get_dimension(nnx,nny,nnz);
	g.get_grid_size(ndx,ndy,ndz);
	g.get_origin(nx0,ny0,nz0);
	if(nnx != ncel_x || nny !=ncel_y || nnz!=ncel_z) return false;
	if(ndx != dx || ndy != dy || ndz != dz) return false;
	if(nx0 != x0 || ny0 != y0 || nz0 != z0) return false;
	return true;
}



GridSpace *GridSpace::divide(GridSpace & g)
{
		GridSpace* newGrid=new GridSpace;
		if(!this->eql(g)) throw GridSpaceNotSame();
		newGrid->set_dimension(ncel_x,ncel_y,ncel_z);
		newGrid->set_grid_size(dx,dy,dz);
		newGrid->set_origin(x0,y0,z0);
		newGrid->initialize(0.0);
		double* d1;
		double* d2;
		d1=newGrid->values();
		d2=g.values();
		for(unsigned long i=0;i<ncel;i++) {
			d1[i]=data[i]/d2[i];
		}
		return newGrid;
}

GridSpace* GridSpace::substract(GridSpace & g) {
	GridSpace* newGrid=new GridSpace;
	if(!this->eql(g)) throw GridSpaceNotSame();
	newGrid->set_dimension(ncel_x,ncel_y,ncel_z);
	newGrid->set_grid_size(dx,dy,dz);
	newGrid->set_origin(x0,y0,z0);
	newGrid->initialize(0.0);
	double* d1;
	double* d2;
	d1=newGrid->values();
	d2=g.values();
	for(unsigned long i=0;i<ncel;i++) {
		d1[i]=data[i]-d2[i];
	}
return newGrid;

}












